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We present a study of a new class of exact solutions having a form of spiral vortices for an isotropic 
O I' two-dimensional Heisenberg ferromagnet using a continuum theory and direct numerical simulations 

D ' of the spin system on a square lattice. We find their features issued from the conservation laws and 

describe their interaction. Reasons behind the formation of the proper spin configurations on a 

■ square lattice are investigated. 

t-H \ 

i 1. I. INTRODUCTION. 

In the last two decades solitons, vortices and other nonlinear excitations in low-dimensional magnets have attracted 
a great interest of researchers. These excitations play an essential role in two-dimensional (2D) magnetism and 
contribute to breaking of the long-range order in 2D magnets. Magnetic vortices are important for the dynamical 
and thermodynamical properties of magnets, for a review see Refsi 2 ^. Predictions of the Belavin-Polyakov theory^ 
for localized structures with a finite energy (instantons) observed much later experimentally^ gave rise to intensive 
"j^ ' investigations of solitons in 2D magnet g&LSiS. 

S Some years ago magnetic vortices have been directly observed in permalloj*i2iii*i2iiii4 and Co magnetic 
nanodotsi^iSiii. High frequency dynamical properties of the vortex state magnetic dots have been probed by Bril- 
louin light scattering of spin waves^ and X-ray imaging techniqusi^. In recent experiments the spin-wave modes 
excited by magnetic field pulses of small litographically define disks with a spin vortex configuration are imaged 
using time-resolved magneto-optic Kerr microscopy^ and phase sensitive Fourier transformation technique 21 . Both 
i ^ i • axially symmetric dynamical modes showing concentric nodes and symmetry breaking azimuthal eigenmodes having 
azimuthal nodes have been observed. An analysis of the time and frequency dependencies^ of the modes demon- 
t-H , strates that for moderate field pulses and large magnetic elements (several tens of microns) the excitation spectrum 
is dominated by magnetostatic modes. However, as noted by authors 2 ^, when the size of the elements is reduced or 
higher modes are excited, the exchange interactions can, in general, no longer be ignored and the dynamic response 
gradually changes from a purely magnetostatic to an exchange-dominated one. One of the aims of the paper is to show 
\ an existence of nonlinear modes both with circular and azimuthal nodes in a 2D isotropic ferromagnet obtained with 

■ an account only exchange interaction. These modes can be observed as spiral-like vortex configurations. Besides a 
' detailed study of the spiral solutions within both XY and Heisenberg models is of interest for the possible applications 

in the physics of liquid crystals, quantum Hall effect and in the study of biological systems featuring self-organized 
spiral structures 2 ^ 2 ^. 

In this paper, we present a study of the spiral vortices in the isotropic 2D ferromagnet using a continuum theory 
and direct numerical simulations of the spin system on a lattice. Our principal concern is to understand in detail 
the structure of the spiral patterns and to find out reasons of their appearance. Real compounds are not ideal 
systems and lattice defects such as impurities, local fields and anisotropy are present in any material sample. The 
effect of nonmagnetic impurities (vacancies) on vortex and vortex-antivortex structures has recently been studied 
for 2D magnetic models with XY symmetry2S*2LS£. These investigations show that an ideal vortex or their pair 
formations are deformed if the vortex centers are near the vacancy. In the continuum limit an account of the 
• i-h , imputities results in logarithmic singularities in the spin field. On the other hand, in the quantum field studies of 
2D 0(N) models enjoying a continuous symmetry these classical configurations with the logarithmic singularities are 
known as " superinstantons"— . To produce these configurations a novel " superinstanton" boundary conditions (SIBC) 
- - ■ were introduced. These consist of Dirichlet conditions on the boundary of the system, and the additional freezing of 
one spin in the center of the sample. It was argued that unlike to standard free, periodic, and Dirichlet boundary 
conditions SIBC do not possess a well defined perturbation expansionist that means that by fixing the spin in the 
center one change a ground state of the system. 

The paper is organized in the following way. In Sec. II the continuum approximation based upon equations of 
nonlinear spin dynamics is presented. We briefly review existing literature and show that the continuum isotropic 
Heisenberg model yields two types of static solutions. A new class of exact solutions of the model, that are local 
minima of the classical energy, is obtained using a special linearization procedure. We find a harmonic function of 
initial dynamical variables obeying the linear Laplace equation. Then, an inverse transformation gives solutions of 
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the initial nonlinear model as functions of the harmonic solutions of the Laplace equation. As a result, new types of 
exact solutions for 2D isotropic ferromagnet are generated where spiral vortices are of special interest. They are likely 
to be relevant for non-perfect systems with defects. Thus, it has been recently shown that vortices are attracted by 
a nonmagnetic impurity^. 

In Sec. Ill we consider numerical simulations of spin configurations predicted by the continuum theory. Firstly we 
investigate spin textures for the planar xy-model with the imposed SIBC. We show that in the sector with topological 
charge q — the ground state is a logarithmic source of the strength a. We find how the strength depends on a turn 
of the fixed spin in regard to spins at the system edges and check the continuum theory prediction for the energy. We 
consider the simplest formation of such kind of logarithmic sources, a pair including two sources of opposite strengths 
a and —a. Then we investigate the case when the structure takes an out-of-plane form. A numerical simulation gives 
the structures which are close to the nodal solutions of Heisenberg model. In the topological sectors q ^ the ground 
state is either a planar logarithmic spiral (xy-model) or a space spiral vortex with an out-of-plane form (Heisenberg 
model). We show that these spin configurations minimizing an energy with imposed SIBC well reproduce features 
predicted by the continuum theory. 

II. ANALYTICAL RESULTS. 

The model to be investigated is the isotropic spin-5 Heisenberg ferromagnet defined by the Hamiltonian 

H = — S ' Jp n SpS n (1) 

p,n 

where S p represents the spin operator at the site p of a 2D square lattice with the nearest-neighbor distance a , and 
Jpn = J3n,p+3 (J > 0) are the nearest neighbor exchange couplings. The non-linear differential equations describing 
the dynamics of the model can be obtained by taking diagonal matrix elements of the equation of motion for the 
raising operator S p = S p + iSp 1 

~i^=[H,S+] (2) 

of the p-th spin in spin-coherent representation |f2) = FJ \0 P ,4> P ) , where < P < n and < <p p < 2ir parametrize 

p 

the spin states on the unit sphere^. For the bilinear Hamiltonian this results in the system for the classical variables 
{9p, <p p } parametrizing the S p spin 



sin t 



,— ^- = -— } J np sin0 n cos 9p cos(ip p - ip n ) + smO p — ^ J np cos0 n , (3) 



= |- ^ J np sin 9 n sin(ip n - <p p ), (4) 

n 

hereinafter, n runs over the nearest neighbors. In the continuum limit we introduce the fields 6(x,y), ip(x,y), which 
are defined in the (x, y)-plane. The equation of motion for static solutions = ^f- — can be obtained by applying 
the continuum approximation to the equation of spin motion on the discrete lattice. This yields 



AO = sin cos (^¥>)' 
V (sm 2 0V(p) = 0. 



(5) 

A remarkable property of these equations is a conformal invariance that allows us to subdivide their static solutions 
into two groups. For the first group of solutions the expression 

80 dtp 80 dip 
dx dy dy dx 
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does not equal to zero. Then we can obtain new types of solutions from already known ones via conformal transforma- 
tions. Indeed, let Oi(x,y),<Pi(x,y) are some particular solutions of Eqs. (JSJ. We may see by direct calculations that 
the fields 8i(ui(x,y),U2(x,y)), (pi(ui(x,y),U2(x,y)) are also the solutions of the same system provided the u± + iu 2 
is an arbitrary analytic function F of the argument x + iy 

ui + iu 2 = F(x + iy). (7) 

Until now, all known solutions belong to the first group. For this case <p{x, y) may be written in the simple form 

(f(x,y) = u 2 (x,y), (8) 

and the another function 8 depends only on the Ui(x,y), i.e. 6i(u\(x, y), u<z(x, y)) — and obeys the simple 

equation of pendulum motion 



9 UlUl (u l ) = -sm[29(u 1 )]. (9) 
The cases of the infinite and finite pendulum motions correspond to the following solutions 



cos#(x, y) = sn 



ui(x,y) ' 



(0 <*<!), (10) 



k 

cosd(x,y) = k sn [ui(x , y) , k] (0 < k < 1), (11) 

where the Jacobi elliptic function of modulus k is used. Eas. (|8ll0lll|l describe a set of quiescent topological defects 
centered at positions z, = + iyi (see Ref^) with 

JL^ / 2i k K \ 

ui(x, y) + iu2(x, iy — Zi) (Ni,Qi€Z). (12) 

i=i V ^ / 

Here, K = K(k) is the complete elliptic integral of the first kind, and Ci are fixed complex parameters. 
For n = 1, Ni = and k = 1 the solution of Eas. 1)811011 

cos8(x 1 y) = taxihui(x 1 y), u\ = Q In \/ x 2 + y 2 , u 2 = Q arctan ^— ^ 

coincides with the Belavin-Polyakov vortex ("baby" soliton)* with the topological charge Q. 

For n = 1, Ni ^ 0, k ^ 1 the solution of Eas. (|10ll2|) represents a A^i-armed logarithmic spiral consisting of 2N\ 
spiral regions separated by the same number of logarithmic spiral walls^. The field tp(x, y) I|8I12|I determined by 
the topological charge Qi does not alter along the logarithmic spiral curves in the (x, y)-plane. Note that the spiral 
vortices in the ferromagnet, involving both 6 and tp variables, have the different mathematical structure in comparison 
with the optical spiral vortices and the spiral vortex solution of the complex Landau-Lifshitz model^ where only ip 
angle is used to build proper configurations. 

We are more interested in the second group of solutions when the expression @ equals to zero. In this case the 
angle ip(x,y) is an arbitrary function of 6(x,y). Then we use the ansatz Vip = f(6)W8 to find solutions of this class 
with V<p\\ V8. After eliminating the V<p from [Eq.JSJl], the equation for 8 being 

sin 6^ + 2/ cos 6> + / 3 sin 2 (9 cos (9 = 0. (13) 
d9 

This is Bernoulli equation in the variable f(9), the general solutions are therefore given by 



f(0) = (W) 
Vc 2 sin 4 8 - sin 2 8 



where c 2 > 1 is an arbitrary parameter. Inserting the ratio Vip — f(8)V8 and Eq. l|14|l into Eq.JSJ) we find the fields 

8(f), ip(f) as 

cos6> = vc c cos a, / 15 s 
(p = arctan (c tana) + ipo. 
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where the field a(x,y) is satisfied the Laplace equation 

Aa = 0. (16) 
We will only consider special case of solution for a(x, y) 



i ( V (x ~ x 0l ) 2 + (y - y 0i ) 2 \ -A / y - y oi \ 

a = at In + } qj arctan I — I , q t e Z. (17) 

i=i V 1 J i=i \x~x 0i/ / 

with the parameters on, g,-, Ri, c. In the above expression, (xw, yoi) and (xoi, yai) are positions of sources and vortices, 
respectively, tpo is an initial value of azimuthal angle tp. We call the parameter a a strength of source that is identical 
to the term used in the hydrodynamic theory. From Eq. l|15fl we see that the parameter c governs out-of-plane spin 
components. In the soliton (|15f) the spins are confined to the vicinity of the ccy-plane with ir /2 — 9 max < 9 < Tr/2+9 max 
unlike the solutions of Eqs. (|8I1UI12I) . The maximal value 9 max is given by 9 max = arcsin Vc 2 — 1/c for n = 1. 
Eas. l|15ll7|l include both novel and well-known solutions considered early by some authors. Below, we list these cases. 

1) n = 1, c= 1 (pure in-plane solutions). Using Ea. i|15fl we find immediately that 9 — ir/2 and ip — q<f> + a\n(r/R) 
written in the polar coordinates r, (f>. For a = 0we restore Kosterlitz-Thouless (KT) vortices, and for q = we have 
"sources", the equation for p being ip = aln(r/R)^^. 

2) n = 1, a ^ 0, q = (solutions with out-of-plane spin components). Eq. i|15fl can be written as 

Vc 2 - 1 / r\ 

cos 9 = cos am — , 

c V RJ ' 



</? = arctan 



(c tan (a In + p (18) 



This agrees with the result obtained in Ref^ ("nodal" solutions). 
A new class of exact solutions 



T 

cos 9 = — cos ( a In — 

R 



1 / 
— cos ( c 



ip = arctan yc tan ya In + 9<^J^ + (19) 

are the two-dimensional spirals^. Figure 1 presents them for different parameters a and q. The former value assigns 
a spiral twist and the latter defines a number of spiral arms. For a = we obtain a vortex with a non-zero out-of 
plane component 

Vc 2 - 1 

cos = cos (q</>) , 

c 

= arctan (c tan ((?</>)) + </?o. (20) 

Unlike the Skyrmion the soliton has a zero topological charge -^(S^) = 0. Since 9 does not depend on the radial 
coordinate r, the solution has no axial symmetry. 

There are several conserved quantities that are important in what follows: the total energy E, the linear momentum 
P, the angular momentum L z and the total number of spin reversals N (see, e.g., Refi^). 

The energy of the soliton given by Eas. (|15|) can be evaluated in continuum approximation resulting in a more 
compact form 




Using Ea. f!7p for the function a we obtain 



E = irJS 2 



E 



In- 



L 



a ; 



qiqj) In 



L 



(22) 
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dij = (xoi — xoj) 2 + (yoi — yoj) 2 is a distance between i-th and j-th vortices, L is a size of the system and ro is a 
cut-off radius where the continuum approximation breaks down. After the formal substitution qiqj — > onctj + qiq^ in 
Ea.l)22[l we recognize the energy of interacting in-plane vortices. We note the absence of terms aiqj. One can see that 
vortices and sources do not interact with each other. The parameters Ri and c do not enter into the expression at all. 
Two spiral vortices of opposite values (a, q) and (—a, — (?) with a pair separation d has the finite energy 

E = 2irJS 2 (a 2 + q 2 ) In ( - 

meaning that these vortices may be bound in pairs. 

Next we derive another conserved quantities related with the spiral vortex. The density of momentum is determined 
by the formula 

* KS - 
P = == Va 

c + v c 2 — 1 cos a 

consisting of the radial 

HS a 

Pr = ■ r^— , (23) 

c+ — lcosa r 

and the azimuthal 

P* = — 7^ Q - (24) 
c + V c — 1 cos a r 

parts. The P r and components are defined by the twist parameter a and the vorticity q, respectively. Substituting 
(J23 E} in the relations 

27T 27T 

f 1 

PrRodcf) — afrS I d<p , — = 2-KahS 

: + v c? — 1 cos (q6 + a In 



o o 
and 

2tt 2tt 



/ 



P^Rodcj) = qhS I d(f> — — — = 2-KqhS 

c + v c — 1 cos ( q<p + a in -^J 



we clarify the physical meaning of the quantities a and q. The first constraint is related with a flow of the momentum 
through the circle of the radius Rq surrounding the vortex core and the last one determines a quantized circulation 
along the circle. 

The total orbital angular momentum L z along the rotation axis through the area irL 2 here reads as 

L z = Shq / rdr I d(f> -, = ShqirL 2 , 

Jo Jo c + Vc 2 — lcos (q4> + aln-^) 

where the density of the angular momentum is defined as 

Shq 

c + Vc 2 — 1 cos (q4> + a In ^) 

The total linear momentum amounts here to J Pd 2 r = 0. 

The conservation of total number of spin reversals (magnon density) 



— cos (a In J — + q<pj J (25) 
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involves the magnon density current 




From Eqs.(J2H| 126(1 . it follows that the presence of additional terms with a non-zero strength of source a produces 
radial components in the densities of the momentum P and the magnon current j. 

To complete our analytical study we discuss stability of the spiral vortices. Given the energy (|21|l we would like 
to consider the effect of perturbation 8ijj, that belongs to the same class as the a(x,y) does, on the soliton structure 
[Eq.O]. This yields 




= - JS 2 j dr ((Va) 2 + 2(Va)(Vty) + (V^) 2 ) = 

= \jS 2 J dr ((Va) 2 + (V^) 2 ) > E[a], 

i.e. the system pays an energy cost and the soliton turns out to be stable against the small perturbations of the field 
a(x,y). 

Among the several questions that should be arisen in the above analysis, an origin of logarithmic sources and a 
range of possible values of a parameter are ones of the most important. Due to the intrinsic interest, an analysis 
of physical reasons behind the formation of spiral vortices is called for. In addition, the continuum theory cannot 
completely describe subtle differences occurring on the lattice at short length scales. Therefore we should concern 
how to organize a numerical process leading to spin configurations that may be compared with the spiral vortices. 
We are aware of the difference between these spin patterns and spiral vortices predicted by the continuum theory. In 
the first case we deal with a ground state of the system under certain constrains and in the other case with stationary 
nonlinear excitations. Our studies will involve lattice model on 2D square lattice. Ultimately, a lattice model is the 
original source of any continuum theoretical description. We will then compare the continuum theory predictions to 
results found in numerical calculations. 



III. NUMERICAL SIMULATIONS 



A. The model. 



To describe in full detail the method of numerical simulations we rewrite the system (|3I4|I in the form convenient 
for an iteration procedure. From Eq.Q we get 



cose 



^ J np sin 9 n sin <p n J = sin ip p I ^ J np sin 8 n cos ip n J 

V n / \ n / 



that yields 



J np sin 9 n sin cp n 
J np sin 6 n sin ip n \ + [Y Jnp sin 6 n cos ip, 



8iny> p = ±-r= (27) 



Jnp sin n cos ip„ 



COS ifip = ±- 



2 / \ 2 

Y Jnp sin 6 n sin tp n I + [J2 J n P sin 9 n cos <p n 



(28) 



7 



and the upper sign must be taken for a ferromagnet (J np > 0). Similar equation for P is obtained from Eq.JSJ) which 
can be written as 

sin 9p I J np cos 9 n ] = cos 8 p J np sin 8 n (cos tp n cos (f p + sin ip n sin ip p ) . (29) 

\ n / n 

Application of i|27l28[l to this equation gives 

y ] Jnp COS On 

(30) 




that after some simplifications yields the expression used in a numerical algorithm 

Jnp cos n 

cos Op = — " (31) 




Together with (|30|l it implies sin(9 p > 0. 

In actual practice, the spin configuration was found by using the original lattice spin fields S n and iteratively 
repointing each along the effective local field due to its neighbors. Scanning linearly through the lattice each site was 
updated in sequence, being reset along the net field due partly to some unchanged neighbors and some that have 
already been repointed. This gives fast convergence than a synchronized global update. The iterations stop if the 
sum 



\ 



N N 

E 

i,j=0 



9 iv 9 

(dnffW - m.< : ) + £ (sin ^ - sin^) (32) 



taken over a quarter of the lattice on the fc-th step is less than tolerance 10 -6 -j- 10~ 10 . We employed the lattice 
coordinates in (|32|l for the notation of site indices. 

The most difficult computational problem in carrying out this program is to find the initial configuration that 
relaxes to a target spin configuration. It is meaningful to impose appropriate boundary conditions too. Obviously 
this a rich problem with a wide choice of options. 

One way is to take the configuration according to continuum formula and assume that each spin has small amplitude 
dynamic deviations from the starting structure. Another approach has been used in a study of a single magnetic 
vacancy centered in vorte:»2&. For numerical calculations a finite core circular system of radius R c is taken. Lattice 
sites are set up surrounding the origin and only those within radius R c are kept. A detailed discussion of starting 
configurations needed for a finding of proper lattice structure and boundary conditions adopted in calculations will be 
given in every case. We note here, they are essentially different for vortex, logarithmic and spiral spin arrangements 
and their pairs. 

Another important point of numerical simulations is a criteria of consistent between the relaxed spin configuration 
and an appropriate continuum solution. We suggest the following scheme for the comparing, (i) The continuum 
theory is not relevant to spins at sites close to the vortex core. Far from the core the relaxed spin angles must well be 
described by the continuum formula. Thus, we have to control this coincidence with a prescribed precision in a region 
where the continuum description works, (ii) The number of independent parameters in a continuum solution must 
be the same as a number of corresponding degrees of freedom controlled in numerical simulations, (iii) A relaxed 
configuration should not lose a symmetry of continuum solution, i.e. it should have a similar dependence on the 
space coordinates (r, </>). (iv) In addition, we confirm the finding of proper solution by analyzing the total energy of 
a relaxed configuration 

S 2 

E = — J np [sill n COS Op COs(tp p — (fin) + COS Op COS n ] 

n,p 



comparing it with a continuum theory prediction. 



8 



B. Logarithmic source 

1. XY-model 

Our current aim is to perform simulations of a logarithmic source in the planar XY model. The static in-plane 
angles satisfy the discrete nonlinear equations 

J2 Jnp Sin (fin 

(33) 




J2 J np COS (fi„ 
Jnp Sin (f n ) + [Yl Jnp COS tp n 



cos ip v = — " = (34) 

1/ \2 / \2 



drawn from Eas. (|27I28|I after the substitution sin0 n = 1. 

A homogeneous arrangement ip p = ip — const is an obvious solution of these equations for any set of the exchange 
couplings J np . An attempt to carry out numerical simulations using a magnetic vacancy with some zero nearest- 
neighbor exchange couplings, a magnetic impurity with another spin and/or different exchange leads to the uniform 
arrangement and does not permit to get a structure similar to a logarithmic solution 

v 

ip = Lp + a\n—. (35) 
K 

A close examination of an in-plane arrangement given by the analytical model shows that there is some bending of spins 
in the center with reference to spin order in the vicinity of the system edges. This nonuniformity of spin distribution 
can gain insight into the numerical process driving the starting spin configuration to desired in-plane structure. A 
reason behind the formation of the nonuniformity may be either a local magnetic field or a local anisotropy. We fix 
our attention on the first case. Indeed, an inclusion into the Hamiltonian of the local field h, acting on a spin at 
position r*o, whose direction relative to fixed spins at boundary is determined by the angle w (see Fig. 2) 

H z = —gnohS J dr 8(f — ro) cos [ui — <p{r)] (36) 

leads to the continuum equation 

A^(f) = ^j^sin[o;- ip(f f )]S(f-f ). (37) 

^•{r) = <Pu(r) + sin [w - tp(f )] In |f - r | , (38) 

where tp u is an arbitrary solution of the Laplace equation Atp u — in two dimensions. A comparison of the result 
with Eq.||33) yields 

a = f^sin[ W -^(r )]. (39) 

We can proceed and estimate the expression using consideration of the mean-field theory. Taking u> = ir/2 and 
suggesting that all the spins are aligned parallel to the one direction except the core spin, pointed fixedly along an 
axis determined both the exchange field zJS of the nearest neighbors and the local filed h, we easily find 

sin [u - ip(r )} = == ==, (40) 

y {zJS) 2 + (g^ohf 



It then follows that 
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where z is the number of nearest neighbors. Therefore, 

a =^—j= = = ■ (41) 

As we can see a ranges from to w 0.636 when h increases from zero to infinity. Most importantly, the local 
magnetic field is the reason of an appearance of the logarithmic source in the system and there is the upper limit for 
a values. 

To check the predictions with the numerical simulation data we consider a square lattice of size (2iV+ 1) x (2A^ + 1) 
shown in Fig. 3 and take J np = J for simplicity. We carry out the iteration process beginning from one of the corners 
of the lattice. Let the local magnetic held directed along j-axis acts on a core spin which has the coordinates (N, N) 
[Fig. 4a]. This spin should be included into the numerical scheme [Eas. (|33l I34|) ] with a little modification 



J] cos Vn + g^h/(JS) 

n 

Y,sin<Pn) + ( Y, cos ¥n + gfih/(JS) 



COS (fNN — — = (42) 

if \2 / \2 



During the iterations a turn of the central spin ipifo) proportional to the applied held h is seen to develop in regard 
to uniform arrangement at the boundaries which hold fixed and are not updated during the iterations [Fig. 4b]. 

It is important to establish that our analytical model reproduces the results obtained by numerical simulations. On 
a lattice the in-plane angles tp n deviate from the formula 1)35(1 and obtain modifications largest near the core spin. 
These angles satisfy a discrete Laplace equation 

^ sin(( y 9„ - ip n+s ) = (43) 

6 

issued from Eq.Q. 

Similar to the analytical prediction, we find that the radial dependence of the arrangement is mostly preserved, 
ip n is determined only by an absolute distance r measured from the core site (N, N) beginning with r ~ 3. Figure 
4c compares the numerical calculations performed for different scan directions with the analytical values given by 
Ea. (|35|l . The points in the inset are fitted by ip = tt/2 + 0.135 In ( 10 r 84 ) ■ The mean-field approximation (|41|) would 
give a = 2/vT77r w 0.154, where the magnetic field is g^h/{JS) = 1. One sees that the agreement is very good. 

It makes sense to explore on a lattice the dependence a[<^(ro)]- Instead of inclusion of local in-plane magnetic field 
the following simple but effective scheme was applied to enforce a desired logarithmic source position. Again a square 
lattice is used with a fixed spin at its center. This core spin with a given tp(r ) would be excluded from updating in 
the numerical routine. In order to find a range of values for the coefficient a we investigate systems of size 21 x 21 and 
41 x 41. Figure 5a summarizes the results found from the relaxed configurations cfij with different starting deviations 
f(ro) by showing a as a function of y(rb)- One sees the dependence is almost linear. The numerical a values for 
N = 21,41 are slightly different, i.e. a weakly depends on the lattice size. The results indicate that relevant a values 
are small, i.e. |a| < 1. More significantly, the analytical model only contains two parameters, a and R, and a relaxed 
spin configuration, obtained numerically, can only depend upon two independent variables <p(fo) and N. 

We then used the numerical code to calculate the energy of the structure 

Jcp = 2 COS fa" ~ + COS fa" ~ ~ E °> ( 44 ) 
np n 

where the first sum runs over all inner sites, the sum in the second term goes over the boundary sites that belongs 
to the structure. The energy was calculated relative to the ground-state energy Eq, an amount of JS 2 per exchange 
bond, for spins aligned with an uniform angle tpo within the xy plane. The energy found as a function of a 2 , taken 
from the fitting, is linear right up to a maximal value a max corresponding to tp{ro) = 7rthat well agrees with the 
continuum model (Fig. 5b). 



2. Pair of logarithmic sources 

It is interesting to confirm the results found above by analysis an assembly of such kind of logarithmic sources. 
According to the continuum model, the simplest formation is a configuration including two sources of strength a and 
—a 
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TABLE I: Data of numerical simulations for pair of logarithmic sources. 



Size E/JS 2 a ipo 



21 x 21 15.9945 0.203 1.562 

51 x 51 15.9984 0.127 1.569 

31 x 31 15.9973 0.160 1.567 

101 x 101 15.9988 0.108 1.5703 

201 x 201 15.9989 0.1005 1.5707 

301 x 301 15.9989 0.1005 1.5707 



a, (g-jgi) +{y-yi) oi (x-x 2 ) + (y-y 2 ) , ,a (x-xi) + {y-yi) , 

^=9 ln P 9 In -J = 0o + - In- -2— -2, (45) 

1 n i z n 2 z {x - X2) + [y - yi) 

where (j)Q = a In [R^/Rx)- The energy of the pair 

E = 2na 2 JS 2 \nd (46) 



is finite and have no dependence on system size. 

We use a square system of size 2N+1 with two sources placed symmetrically near its center which has the coordinates 
(TV, N) . We found in the numerical studies that only distance between the sources and their mutual orientation on a 
lattice affect the values of energy and a. We present here our results obtained from simulations when the sources are 
placed at positions (N — d/2, N) and (N + d/2, N), where d is a distance between them. Alternatively, if the sources 
were placed at different positions away from the system center a boundary energy that changes significantly with the 
pair positions would result. To avoid this complication, it is much simpler to fix the source positions at the system 
center. 

We investigate how the strength of source a could depend on the pair distance d and on a difference between the core 
spin angles Sp — p(r2o) — p(fio) an d how the analytical predictions for the energy (E ~ a 2 , E ^ hid) are modified on 
a lattice. The calculation required a larger system than for the previous case in order to produce stable configuration 
in a finite scale. We checked the finite size effects by simulating 21 x 21 to 1001 x 1001 square lattices. In addition, in 
the calculations presented here we control a restoring of homogeneous arrangement at system edges. The analytical 
solution (|45|l involves two independent parameters, i.e. the strength a and the distance d. In numerical simulations 
a number of independent quantities involved in the calculation is the same. The first parameter is determined by the 
difference Sip and the second one is given explicitly via the source positions. 

It is of importance to determine an optimal size of the system to avoid impractically long simulations for large 
enough lattice. In the Table I some of the data, answering the purpose, is collected. The binding energies, the 
strength of source a and a relaxed background arrangement ipo are summarized there. As a starting configuration 
we take two spins at the distance d — 4 turned almost oppositely to another magnetic moments. We run then 
iteration procedure to obtain equilibrium configuration. Fitting to the known solution l|45|l allows the determination 
of the needful quantities. One see that the background arrangement po does not approach the exact value 7r/2 when 
decreasing the size of the lattice (L < 101). Size effects are noticed for the strength of source a which becomes greater 
for small lattices. Increasing the L from 201 till 301 has an effect neither on a nor on po. As a rule of thumb, we 
hold the finite-size effects are negligible when L > 50d. 

To compare the analytical expression l|45l) with a target in-plane arrangement (Fig. 6) we choose a scan along a 
path in the (0, 1) direction of the lattice, beginning with the point of one of the source. Similar results could be 
obtained along other scan directions, but with a worse consistence with the analytical predictions, since the region, 
where the continuum model is expected to be valid, is not isotropic. A fit of data points according to Eas. l|45l I4rj|) 
provides estimates for the parameter a and the energy E. 

The dependence a on Sip was examined for several d values. We found that the observed dependence is distinguished 
for small (d = 2 4) and large (d > 6) pair distances. 

For sufficiently small distances (d < 4) the pair presents an unit formation and cannot be treated adequately by 
the continuum model. A typical result for the energy obtained for a square lattice of size L — 141 is shown in Fig. 
7. The parameter a found as a function of Sip is shown in Fig. 7a. One sees that the dependence exhibits a clear 
periodical behavior. Due to the fact, the dependence E(a 2 ) calculated with the same data and plotted in Fig. 7b 
turns out to be nonlinear and many-valued. 
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TABLE II: Data of numerical simulations for nodal states. 
Sip = ip c — ipb 59 — C — ^ A(, c E/JS 2 c ipo R a a (plane) E (plane) 



1.5 0.5 1.509 1.513 1.140 -f 21.54 -0.317 -0.318 1.513 

1.5 -0.5 1.509 1.513 1.140 -f 21.54 -0.317 -0.318 1.513 

0.5 1.0 1.077 0.744 3.399 § 29.62 -0.227 -0.228 0.714 

2.5 1.5 1.627 1.758 23.58 f 21.69 -0.340 -0.340 1.757 

1.509 1.509 1.513 oo - 21.54 -0.317 -0.318 1.513 



By contrast, for the large distances (d > 10) the agreement between the analytical model and the numerical 
simulations (L = 701) is good enough. In Fig. 8a we show a as a function of Sip. The dependence is almost linear 
and this behavior is similarly observed for one logarithmic source. The energy of the pair as a function of a 2 also 
supports the agreement, E is directly proportional to a 2 until Sip ~ it (Fig. 8b). 



3. Nodal solutions of Heisenberg model. 

The patterns which we have so far studied have been confined in the xy-plane. A particularly interesting and 
complex case is when the structure takes an out-of-plane form. Here we report numerical simulations of nodal states. 
The lattice is taken as in Fig. 3, however, in the initial configuration the pinned spin in the center has an out-of- 
planc component. The numerical procedure recurs those used for one logarithmic source and generates both a set 
(cos ipij , sin ifij ) and (cos 0y , sin 6ij ) . This allows us to extract from these simulations parameters involved in analytical 
expressions such as those derived in Sec. I. 

Given a set of angles found numerically, we should summarize the data by fitting it to the model 118|) that depends 
on the adjustable parameters ipo, c, a and i?, and the fit supplies the appropriate coefficients. To carry out the 
comparison we follow two steps. From Eas. (|18|) we obtain 



\/c 2 — 1 cos ipo cos ipij + \/ c 2 — 1 sin ipo sin tpij = tan j . (47) 

Fitting Ea. H47|l to the resulting spin structure by the least-square method, we obtain c and ipo- The rest parameters 
a and R are then adjusted to data for cos 9 . We estimate 



ipij = tan 1 c tan (^a In ^ 



- fo 



with found a, R, c and ipo and compare ip^ with numerical data points (Fig. 9). One sees, the agreement is rather 
well. 

A set of numerical simulations was performed using different initial conditions for boundary Sf, = 
(sin 0;, cos </3f,, sin #f, sin cos 0;,) spins and the central pinned spin S c — (sin# c cos<p c ,sin# c sinc/j c ,cos# c ). We also 
used the expressions <|18H as fitting formulae. These fits give us values for the strength a and the energy E listed in 

Table II. We found a significant signature: the a value is determined only by the angle At c = cos -1 (^SbS c ^j for any 

direction S c . To check this assertion we repeat calculations of Sec.(III.B.l) for the pure in-plane case with the starting 
deviations ip(ro) = A(, c . A fit of data points according to Eas. l|35l44|l gives the estimates for strength of source a 
and the energy E (last two columns in the Table II) which are close to the out-of-plane values. In addition, these 
simulations confirm the analytical results for the energy: E ~ a 2 and E does not depend on c. 

From these features we could conclude that a pinned spin in the center is a reason of appearance of logarithmic 
solutions both in XY and Heisenberg models. 



C. Spiral vortex 



In the present section we continue with simulations of a spiral vortex, and begin with XY case. It is meaningful 
to investigate whether the approach described earlier keeps its validity for simulations of an ideal vortex. To perform 
numerical simulations we consider a square lattice of even size 2N x 2N shown in Fig. 10a and place the vortex core 
in the center of the dual lattice. Then, we start with an initial configuration with q =/= and impose free boundary 
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TABLE III: Data of numerical simulations for spiral vortex. 



(go, yo) (E - E vortex )/JS 2 a \a\ (energy) 

(N,N) 0.114 -0.084 0.0828 

(N + 5, N) 0.554 -0.179 0.183 

(N + W,N) 0.562 -0.182 0.184 

(AT +15, AT) 0.568 -0.189 0.185 



condition. After 5000 iterations we reach a relaxed configuration shown in Fig. 10b. On a lattice, the in-plane vortex 
angles ipij lose the perfect circular symmetry of this formula, and obtain modifications largest near the vortex core 
with the coordinates (io, Jo)- To check a validity of the solutions 

*-*o . 3~3o ,,m 

cos <fij = — = , sin ipij = — = (48) 



'(i -*o) 2 + io) 2 \f(i-io) 2 + (j-3o) 2 

we control a fulfillment of a discrete lattice nonlinear Laplace equation 

sin {<Pp - <Pn) = 0, 

n 

where n runs over the nearest neighbors of the p-th site, with an accuracy of order 10~ 2 in the center and 10 -5 at 
the outskirts. 

As a next step, we perform numerical simulations for a spiral vortex determined by the most common expression 
for the harmonic function 



a(x, y) = q arctan + a In V^E^EEEM 

t — nr. _ hi 



X X q fi, 

Here we briefly review how this procedure has been organized. Assuming the vortex is placed at some position centered 
in a plaquette and a logarithmic source position coincides with one of the lattice sites the parameters x ,y (so, 2/0 ) 
are chosen in the dual (direct) lattice. Firstly, we consider an initial vortex configuration, where one of the spins 
nearest to the vortex is considered as a site of logarithmic source. The in-plane angle of this core site is not changed 
by the iteration procedure, and a fixed boundary condition holds. Farther from the core the relaxed spin angles must 
well be described by the continuum formula ip(x,y) = q4> + a\n{r/R), then we expect the spin configuration 



j-N +1/2 y/(i - x Q ) 2 + (j - yo) 2 

tfiij = q arctan — — + a In ■ 



i- N +1/2 R 
would be seen to develop. 

We investigate system of size 101 x 101. The starting deviation ip{fo) of the core spin in the uniform background 
results in the relaxed configuration with a and the energy E (Fig. 11a). Most importantly, we observe a logarithmic 
dependence of the in-plane angles ip^ for scans (used further in the fitting to determine a values) along paths in definite 
directions (Fig. lib) while scans along another directions shows the opposite feature (Fig. 11c). This procedure has 
been repeated at several different positions of the fixed spin (xq, yo) and the results are shown in Table III, where in 
the last column we show the a values derived from the calculated energy E. 

This is done by using the expression E = E\ og + E vort ,where 

4 R = nJS 2 a 2 In — , E vort = nJS 2 q 2 In — , 
a a 

for the spiral vortex energy in the continuum approximation which is suggested to be valid. The value for E vort / JS 2 = 
16.58 is taken from the simulations of the ideal vortex. With this general expression, where E vort /E\ og = q 2 /a 2 , it 
is easy to find a which will be given by a (energy) = \/{E — E vort ) /E vort . The following features are evident from 
Table III. 

1) The farther the fixed spin from a vortex center the closer an additional energy E — E vort (second column) to 
value found for a fixed spin in the uniform background. The region of continuum theoretical description decreases 
that is in agreement with the simulations for two opposite logarithmic sources already reported in Sec. III. B. 2. 



13 



TABLE IV: Data of numerical simulations for space vortex. 



c 




E 


1.0 
1.5 
3.5 




2.107 
2.107 
2.107 


TABLE V: Space vortex energy as a function of the inner radius Ri. 


Ri 


E (lattice) 


E = 7T JS 2 In §i 


10.5 
20.5 
30.5 
40.5 
50.5 


7.05 
4.93 
3.68 
2.80 
2.11 


7.10 
5.00 
3.75 
2.86 
2.16 



2) The strength a depends both on f(ro) value and a position of the fixed spin. The closer latter to a vortex center 
the less a at fixed <p(ro) value. Far from the vortex source the parameter a becomes exactly equal to value found for 
one logarithmic source with the same lattice size and starting deviation <f{fo). 

3) The prediction of continuum theory E ~ (a 2 + g 2 ) is full reproduced by simulations for any logarithmic source 
and vortex positions. The lowest energy value is obtained for the smallest distance between the vortex center and the 
fixed spin. 

D. Space spiral vortex 

We begin with the case a = 0. To obtain this kind of spin structures we solved the equations 1)271 1281 l3*T|l by first 
setting the angles to their continuum values l|20l) and then iteratively setting each spin components to point along the 
direction of the effective field due to its neighbors. However, an attempt to obtain the space vortex by this way on 
the full square or disk fails. For these systems the iteration procedure either converges to a Skyrmion structure (free 
boundary conditions) or does not converge at all when the boundary spins held fixed. The reason is the iteration 
procedure relaxes to a minimal energy state. Since the feature of a starting configuration employed in numerical 
calculations is a non-zero angular momentum, this configuration may evolve either into the Skyrmion-like or into the 
space vortex-like structures. However, the former has a gain in energy in comparison with the last one. To avoid this 
difficulty we consider a simple way in which the Skyrmion structure loses the advantage. Noting that an essential 
contribution to the vortex energy comes from the spins within a small-radius core we takes the computional region in 
the form of ring with the inner radius R\ and the external radius R2 . This results in small differences of the discrete 
solution from the continuum result 12UI) . 

We found the energy for the different parameters c as shown in Table IV for the ring of size Ri — 50.5 and 
i?2 = 105.5 with free boundary conditions. 

In full agreement with the continuum theory these energies have no dependence on c values. In Table V we listed 
the space vortex energy as a function of the inner radius R\ at fixed R2 — 100.5. We see that agreement between the 
numerical and the continuum theory result E = ttJS 2 In is nice for the whole range of radius 10.5 < Ri < 50.5. 
We may conclude that in the system without a finite-radius core the space vortex will have the lowest energy among 
solutions with a non-zero angular momentum. This assertion is supported by direct analytical consideration^. 

Now we turn to the space spiral vortex with a 0. To obtain this configuration spins belonging to both boundaries 
of the ring are held fixed. The constants <po(Ri) and 0o(i?2) should be taken different ("twisted" boundary conditions). 
Figure 12 presents the results of the modeling for R\ — 50.5, R2 = 105.5 and 4>o{R2) — <Po(Ri) + 7r/ '4 resulting in a 
spiral vortex structure with the energy E = 2.37 JS 2 . 

Recently, the statics and dynamics of flat circular magnetic nanostructures with an in-plane magnetic vortex 
configuration has been investigated within the framework of the Landau-Lifshitz-Gilbert equation putting particular 
emphasis on the polarization of the vortex center and on the in-plane vorticity^S. Studying fast switching process 



14 



induced by out-of-plane field pulses, the authors was no longer dealing with a vortex state, but rather with a spiral 
(see figure 12 in 40 ). They found also that in nanorings with an inner radius R\ and an outer radius R2 the stability 
of the vortex state is enhanced, and concerning the switching of the vorticity, the nanorings have similar properties 
as circular ones, i.e. with R\ — 0. 

In summary, we have studied a new class of spiral vortex-type solutions in a 2D Heisenberg ferromagnet and 
performed numerical simulations for various spiral vortex configurations using fixed twisted boundary conditions and 
pinned core spins (" superinstanton" boundary conditions). These simulations show a reasonable agreement with 
the continuum- approximation results. Based on the investigation we may identify among the nonlinear excitations 
the modes with circular nodes ("nodal" solutions^) and modes with azimuthal nodes of magnetization M z (space 
spiral vortex) that resembles the classification of magnetostatic modes excited by a magnetic-field pulses and observed 
recently in micron-sized ferromagnetic disks. Incidentally, only axially symmetric magnetostatic modes appeared if the 
tipping pulse is uniform over the disk and all geometries are perfectly axially symmetric. Symmetry breaking modes, 
instead, required, e.g., a non- uniform tipping pulse having a sizable gradient in the plane of the vortex or a deviation of 
the sample from a perfect cylindrical shaped. However, the frequency of nonaxially symmetric magnetostatic modes 
has a negative dispersion, i.e. it decreases with a growth of a number of azimuthal nodal lines. Unlike this, for the 
nonlinear excitations, which are of exchange origin, a number of these lines coincides with a number of spiral arms 
(or with the vortex topological charge q) and increases the energy. These stationary nonlinear modes must be taken 
into account for yielding a better understanding, e.g., the fast magnetic switching properties in magnetic memory 
materials. 
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FIG. 1: Spiral vortices with a — 2.0 [q — 1 (a), q — 2 (b)] and with a — 0.3 [q = 1 (c), q — 2 (d)]. Density images of 
the amplitude and the phase of the magnetization are given for the last set. The map (e) has one node curve dividing the 
bright-dark regions of equal amplitude but opposite phase. The map (f) consists of four regions, oscillating in pairs in phase. 
Note a similarity of topological small-a spiral vortices with dynamical magnetostatic nonaxially symmetric modes (see Fig. 1 
in Ref. [23]) and Figs. 3-4 in Ref. [21]. 



FIG. 2: Core spin (white arrow) points along an axis determined jointly by the exchange field zJS of the nearest neighbors 
and the local filed h. 



FIG. 4: Starting configuration used to get a logarithmic source (a) and equilibrium configuration obtained in the iteration 
procedure (b). Core spin is placed at position (10, 10). The final equilibrium state reached after 50000 iterations with an 
accuracy 10~ 15 . A comparison between the numerical simulation (points) and the analytical model (lines). To exclude the 
lattice artifacts due to fixed boundary conditions we use the points with 1 < lnr < 2 for the fitting shown in the inset (c). 



FIG. 5: The a value as a function of <p{ro) (a); the energy E as a function of a 2 (b). Note that the maximal a values found for 
different lattice sizes, a ma x = 0.650 (N = 21) and a max = 0.615 (N = 41), are close to the estimation 0.636 of the mean-field 
theory. 



FIG. 6: Starting configuration used to get a pair of logarithmic sources (a) and relaxed configuration obtained in the iteration 
procedure with an accuracy ~ 10 -7 (b). In the starting configuration the in-plane angles of the selected sites are set equal to 
7r/2 + Sip/2 and n/2 — 5<p/2 and they are hold fixed during the iteration scheme [Eqs. H33l34ll 1. The spins at the edges of the 
system have been included into the iteration procedure too and they have no constraint from outside, that is a free boundary 
condition holds. For all another spins, the in-plane angles are supposed to be it/2. 



FIG. 7: Parameter a as a function of Sip (a) and energy £ as a function of a 2 (b) for a small-distance pair (d = 4). The pair 
presents an unit formation and cannot be treated adequately by the continuum model. 



FIG. 8: Dependencies a(Sip) (a) and E(a 2 ) (b) for a large-distance pair (d — 10). The agreement between the analytical model 
and the numerical simulation is good enough. 



FIG. 9: Comparison of cos<£ obtained from numerical simulation (black points) to the analytical expression cosi^ = 
cos {arctan [1.140 tan (—0.307 lnr + 2.544)] — n/2} given by the continuum theory (solid line). Core spin angles are <p(ro) — 1.5 
and 9(r ) = 0.5. 



FIG. 10: Starting configuration for a simulation of vortex structures on a lattice of size (2N, 2N) (a). The center with 
coordinates (N — 1/2, N — 1/2) is denoted by black circle. Relaxed configuration obtained after 5000 iterations (b). The 
estimation of energy E/JS 2 for system of size 100 x 100 by using numerical codes ipij [Eg. 1441 1 is 16.58 vs 17.72 given by the 
continuum theory. 



FIG. 11: Relaxed spiral-vortex configuration with the energy E/JS 2 = 16.6994. Pinned (N, 7V)-spin is denoted by the dotted 
circle (a). Dependencies of the in-plane angles <pij along the ( — 1,-1) (b) and (1,-1) (c) directions with logarithmic and 
non-logarithmic behavior, respectively. 



FIG. 12: Space spiral vortex (c=2): arrangement in the plane xy (a) and the cos# profile(b). 



FIG. 3: Coordinates for square lattice. The core spin is denoted by white-black circle. Solid circles indicate the inner 

sites involved into iteration procedure. Open circles are the boundary sites. 
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